include	<error.h>
include	<imio.h>
include	<imhdr.h>

task datacube

procedure datacube()

define MAXLENS 2000

pointer outim,inptr,image,imout,im,outptr
int imlist
int outlist
int x,y
int nchan,nesp

char ldispdir[SZ_FNAME]
char lconf[SZ_FNAME]
char lconfpath[SZ_FNAME]
pointer fd
int  xl[MAXLENS]
int  yl[MAXLENS]
int  nl
int  xsize, ysize
int stat
int i,j

pointer open (), fscan (), close ()
pointer ahivi()
int	imtopenp(),imtgetim()
pointer	immap(), imgs2r(), imps3r()


begin

	call salloc (image, SZ_FNAME, TY_CHAR)
	call salloc (imout, SZ_FNAME, TY_CHAR)

	imlist = imtopenp ("image")
	outlist = imtopenp ("imout")
	call clgstr ("ldispdir", ldispdir, SZ_FNAME)
	call clgstr ("lconf",lconf, SZ_FNAME)


	call strcpy(ldispdir,lconfpath,SZ_FNAME)
	call strcat("/",lconfpath,SZ_FNAME)
	call strcat(lconf,lconfpath,SZ_FNAME)
	fd = open (lconfpath, READ_ONLY, TEXT_FILE)
	nl = 0
	while ( fscan (fd) != EOF ) {
        call gargi (xl[nl])
        call gargi (yl[nl])
		nl = nl + 1
	}
	stat = close (fd)

	if (nl == 0) {
		call error (1, "Can't find any lens positions in the config file.\n")
	}

	xsize = ahivi(xl,nl) + 1
	ysize = ahivi(yl,nl) + 1

	call printf ("Found %d lens positions in the file '%s' \n")
		call pargi (nl)
		call pargstr (lconf)

	call printf ("IFU dimensions: %d x %d \n")
		call pargi (xsize)
		call pargi (ysize)

	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF && imtgetim (outlist, Memc[imout], SZ_FNAME) != EOF ) {
	    iferr {
			im    = immap (Memc[image], READ_ONLY, 0)
			outim = immap (Memc[imout], NEW_COPY, im)

			nchan = IM_LEN(im, 1)
			nesp  = IM_LEN(im, 2)

			call printf ("Found %d spectra each with %d channels \n")
				call pargi (nesp)
				call pargi (nchan)

			call imputi (outim, "i_pixtype", TY_REAL)
			call imputi (outim, "i_naxis", 3)
			call imputi (outim, "i_naxis1", xsize)
			call imputi (outim, "i_naxis2", ysize)
			call imputi (outim, "i_naxis3", nchan)
			
			inptr = imgs2r(im,1,nchan,1,nesp)
			outptr = imps3r(outim,1,xsize,1,ysize,1,nchan)
			
#			call amovr(Memr[inptr],Memr[outptr],nchan*nesp)
			
			do i=0, nesp-1 {
				x = xl[i]
				y = yl[i]
				do j=0, nchan-1 {
					Memr[outptr+x+y*xsize+j*xsize*ysize] = Memr[inptr+j+i*nchan]
				}
			}
			
			call imunmap (outim)
			call imunmap (im)

		} then {
			call erract (EA_WARN)
			next
	    }
		
		
	}

end

